%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Program to accompany "Using VARS to Identify Models of Fiscal Policy: A
%Comment on Perotti," NBER Macroeconomics Annual 2007, by Ricardo Reis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This file generates the figures in the comment. See the accompanying
%Mathematica file for the proofs of Propositions 1 and 2.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear; clc; close all

%Setting structural parameters.
global beta theta csi alpha delta sigma GY taoK s YK GK CK Xi
beta = 0.99; theta = 1; csi = 4; alpha = 0.34; delta = 0.025;
GY = 0.21; s = 0.12;  taoK = 0.54;

%Defining a set of useful auxiliary parameters
sigma = s/GY;
YK = (inv(beta*(1-taoK))-1+delta)/(alpha*(1-s));
GK = GY*YK;
CK = YK - GK - delta;
Xi = YK*(1 - s)/(YK*(1 - s) + 1 - delta);

%Picking length of impulse repsonses
global T; T = 16;

%Load IRFs estimated using data from 1947
global t_irf g_irf t_sd g_sd
g_irf = [0.057519693	0.066941269	0.074274742	0.082844246	0.080083129	0.075823458	0.069192346	0.061815999	0.054276898	0.046760703	0.039875228	0.033490277	0.027642511	0.022421336	0.017793701	0.013753809	0.010291268	0.007375029	0.004966145	0.003015696	0.00147116	0.000279171	-0.000613235	-0.001255155];
g_u	= [0.059733895	0.066501076	0.072908883	0.080185481	0.074852571	0.068911764	0.060330434	0.051964461	0.043701462	0.035377579	0.028176793	0.020841203	0.014719181	0.009403521	0.004674016	0.000423818	-0.002377478	-0.005533474	-0.007294727	-0.009095346	-0.010511545	-0.01073601	-0.011645157	-0.012663188];
g_l	= [0.066038653	0.079062712	0.090124791	0.100991065	0.100656384	0.097351985	0.092280454	0.085645574	0.078246548	0.070728741	0.062891901	0.055670601	0.049666127	0.043825737	0.038939753	0.034351798	0.029866967	0.027152179	0.024153048	0.020973205	0.018335416	0.015903962	0.013970346	0.012542939];
y_irf = [0.695010851	0.778111428	0.815787513	0.905267032	1.030654086	1.109637457	1.15192838	1.173581757	1.144856401	1.06450403	0.952052294	0.821460054	0.683138935	0.54975957	0.428740535	0.323809263	0.23646448	0.166250916	0.111868319	0.07149858	0.043176006	0.024988035	0.015052882	0.011590769];
y_l = [0.517386677	0.462491636	0.403520548	0.432487235	0.530204009	0.582402222	0.621832755	0.648094959	0.613234551	0.505189935	0.351399846	0.211248935	0.063174743	-0.055817576	-0.191243196	-0.29844099	-0.394438224	-0.433102147	-0.469719392	-0.481238452	-0.509698	-0.527272574	-0.550966478	-0.512855468];
y_u = [0.977501451	1.222647686	1.403081687	1.563203623	1.723800534	1.856025234	1.962438444	2.057724072	1.990593112	1.952464889	1.833491421	1.701778445	1.599861867	1.452317809	1.310843426	1.146269054	1.002693595	0.898083031	0.801875651	0.734136628	0.683058992	0.653902431	0.643124192	0.624012543];
c_irf = [0.001187467	0.003035812	0.002306913	0.004208314	0.005201792	0.005976273	0.006947678	0.007484585	0.007780427	0.007859715	0.007724498	0.007434617	0.007038523	0.006585276	0.006106667	0.005629649	0.00517321	0.004748208	0.004360504	0.004012774	0.003704876	0.003434801	0.003199313	0.002994319];
c_l = [-0.000105504	0.001381621	-2.6986E-05	0.00189154	0.002509706	0.003040042	0.003842694	0.004349933	0.004691459	0.004585193	0.00434291	0.004021164	0.003521339	0.002868381	0.002303581	0.001991279	0.001230918	0.000838332	0.000688896	0.000373752	-0.000148913	-0.000476007	-0.000652521	-0.000845487];
c_u = [0.002760994	0.005480954	0.005335127	0.007555709	0.009227011	0.01042308	0.011856054	0.012695833	0.013301608	0.013612622	0.013691999	0.013669729	0.01339738	0.012999368	0.012395937	0.011560229	0.01105693	0.010581774	0.010079859	0.009637582	0.009106204	0.008714122	0.008197826	0.007847243];
t_irf =[0.325340011	0.502082022	0.331463856	0.523682923	0.467366941	0.461678102	0.460014324	0.413694346	0.376588631	0.328301542	0.279382747	0.231245345	0.183663887	0.141508554	0.104077337	0.071802687	0.045247219	0.024012249	0.007760722	-0.003904591	-0.011484403	-0.015465472	-0.016373754	-0.014725955];
t_l	= [0.208810317	0.351035402	0.125057541	0.302557532	0.211666151	0.204002387	0.157421695	0.124439344	0.071226194	0.009561347	-0.038858494	-0.092622467	-0.13310375	-0.162362458	-0.19492327	-0.226963571	-0.253327215	-0.279162571	-0.2891734	-0.295566517	-0.286889048	-0.281625379	-0.260995326	-0.239330657];
t_u	= [0.500578701	0.754059827	0.613966101	0.848352574	0.800617569	0.817043934	0.818076406	0.799222185	0.772419777	0.715302071	0.658879129	0.601731693	0.543191257	0.502324929	0.478736641	0.433757932	0.403415565	0.37935171	0.348618061	0.334355835	0.312722488	0.291471886	0.276397651	0.272162977];
%Reshape all to be of length T
y_irf=y_irf(1:T); c_irf=c_irf(1:T); t_irf=t_irf(1:T); g_irf=g_irf(1:T);
y_l=y_l(1:T); y_u=y_u(1:T); c_l=c_l(1:T); c_u=c_u(1:T);
t_l=t_l(1:T); t_u=t_u(1:T); g_l=g_l(1:T); g_u=g_u(1:T);

%Calculate parameters that give best fit between model and estimates for
%the impulse responses of government spending and taxes.
options = optimset('MaxFunEvals',100000000);
starting_guess = [0.5 0.1 0.1 0.1 0.1];
[est,loss] = fminunc(@RPfigures_c,starting_guess,options)
disp('Parameters found are: ')
disp('     rho     gammaG    gammaT     lambdaG  lambdaT: '); disp(est);
disp('And the SSR is:'); disp(loss);
rho = est(1); gammaG = est(2); gammaT = est(3); lambdaG = est(4); lambdaT = est(5);

%Given these parametrs, now solve for policy functions for k and c and for
% their ARMA structure:
[eta H Ac Bc] = RPfigures_a(est);
disp('Check that roots eta and rho are smaller than one for stable system:')
disp([eta rho])
%Solve for policy functions and ARMA for g, t, n and y
Ag = gammaG*Ac;
Bg = gammaG*Bc + (1+lambdaG);
At = (gammaT + lambdaT*gammaG)*Ac;
Bt = (gammaT + lambdaT*gammaG)*Bc + lambdaT*(1+lambdaG);
An = (alpha - Ac*theta - At)/(alpha+csi);
Bn = (-Bt - theta*Bc)/(alpha+csi);
Ay = (1-s)*alpha + (1-s)*(1-alpha)*An + s*Ag;
By = (1-s)*(1-alpha)*Bn + s*Bg;
%Solve for ARMA representations
muc = Bc; nuc = Ac*H - Bc*eta;
mug = Bg; nug = Ag*H - Bg*eta;
mut = Bt; nut = At*H - Bt*eta;
mun = Bn; nun = An*H - Bn*eta;
muy = By; nuy = Ay*H - By*eta;
%Compute predicted IRFs
irfg = RPfigures_b([eta rho mug nug],T);
irft = RPfigures_b([eta rho mut nut],T);
irfc = RPfigures_b([eta rho muc nuc],T);
irfn = RPfigures_b([eta rho mun nun],T);
irfy = RPfigures_b([eta rho muy nuy],T);

%Plot IRFs for t and g
figure
subplot(1,2,1), plot(1:T,t_irf,'-b',1:T,irft,'-or',...
    1:T,t_l,'b:',1:T,t_u,'b:',0:T,zeros(1,T+1),'k','LineWidth',2)
title('Taxes'); legend('Data','Model')
subplot(1,2,2), plot(1:T,g_irf,'-b',1:T,irfg,'-or',...
    1:T,g_l,'b:',1:T,g_u,'b:',0:T,zeros(1,T+1),'k','LineWidth',2)
title('Government spending'); legend('Data','Model')

%Plot IRFs for c and y
% I re-scale them to make the graph more readable.
irfc = irfc.*(sum(abs(c_irf))/sum(abs(irfc)));
irfy = irfy.*(sum(abs(y_irf))/sum(abs(irfy)));
figure
subplot(1,2,1), plot(1:T,y_irf,'-b',1:T,irfy,'-or',...
    1:T,y_l,'b:',1:T,y_u,'b:',0:T,zeros(1,T+1),'k','LineWidth',2)
title('Income'); legend('Data','Model')
subplot(1,2,2), plot(1:T,c_irf,'-b',1:T,irfc,'-or',...
    1:T,c_l,'b:',1:T,c_u,'b:',0:T,zeros(1,T+1),'k','LineWidth',2)
title('Consumption'); legend('Data','Model')
